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Replica symmetry breaking in an adiabatic spin-glass model of adaptive evolution 
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We study evolutionary canalization using a spin-glass model with replica theory, where spins 
and their interactions are dynamic variables whose configurations correspond to phenotypes and 
genotypes, respectively. The spins are updated under temperature Ts, and the genotypes evolve 
under temperature Tj, according to the evolutionary fitness. It is found that adaptation occurs at 
Ts < Tg^^, and a replica symmetric phase emerges at yj*^^ < Ts < Tg^. The replica symmetric 
phase implies canalization, and replica symmetry breaking at lower temperatures indicates loss of 
robustness. 
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Biological evolution occurs through changes in geno- 
types and phenotypes over generations, driven by ran- 
dom genetic variance and natural selection. This process 
preferentially selects genotypes that produce a phenotype 
that affords high evolutionary fitness [H, [I]- Thus, phe- 
notypes, such as protein expression levels or the func- 
tional structures of proteins, are the result of dynamic 
processes governed by the genes. However, such processes 
generally involve stochasticity due to thermal noise, and 
thus phenotypes of isogenic individuals are not necessar- 
ily identical [3|-[5|- Indeed, such phenotypic fluctuations 
and the possible role of noise have been extensively in- 
jgsti^ated both experimentally 0] and theoretically 

For a phenotype to conserve its function, however, it 
must be robust to this noise, at least to some degree. In- 
deed, the dynamic adaptation process that shapes pheno- 
types exhibits global and smooth attraction, as observed 
in the folding dynamics of proteins [ill, [13 , R-NA [l^ , 
protein expression dynamics governed by gene regulatory 
networks [l4| , developmental dynamics |15| , and so forth. 
Besides this robustness to noise, the adapted phenotype 
should be robust to genetic change to acquire evolution- 
ary stability. The possible relationship between these two 
types of robustness, as well as the positive role of noise, 
has recently been investigated theoretically [sl-fiol [Tsl. [l6j . 
The study found a transition toward robustness in the 
dynamic process with respect to the noise level (tem- 
perature), where the energy landscape for the dynamics 
changes from being rugged to having a funnel-like struc- 
ture. 

Considering the above change in the dynamical pro- 
cess, one may expect that loss of robustness could be 
viewed as a transition to the spin-glass phase in statis- 
tical physics. Thus far, however, no analytic theory to 
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support this view has been provided, and, from a theoret- 
ical standpoint, little is understood of this transition in 
the evolution of robustness against noise (temperature) . 

Here we introduce a simple statistical-mechanics model 
of adaptive evolution to explain the dynamical process 
that shapes phenotypes. We use an adiabatic two- 
temperature spin-glass model in which the spin config- 
uration and the interaction matrix correspond to the 
phenotype and genotype, respectively. The genotype 
evolves to increase fitness which is defined by the spin 
configuration. With an analysis based on replica theory, 
we demonstrate the emergence of a replica-symmetry- 
breaking transition as the temperature decreases, and 
show that the transition corresponds to a loss of ro- 
bustness in the phenotype. Adaptive evolution of ro- 
bustness is shown to occur only in the replica symmetry 
phase, where the Hamiltonian for global attraction to 
the adapted phenotype is represented in terms of frustra- 
tion. We also discuss the significance of replica symmetry 
breaking on phenotypic robustness. 

Let us consider a simple spin model in which the phe- 
notype and genotype are represented by configurations 
of spin variables S = {Si} and the interaction matrix 
elements J = {Jij}, respectively, with i,j = I,-- - ,N. 
Each spin variable Si can take one of two values, ±1. 
The interactions are fully connected between two spins. 
Both the spins S and interactions J arc treated as dy- 
namical variables, but the time scale associated with J is 
much slower so that the interactions arc relatively fixed 
during the time evolution of the spins. Thus, the equi- 
librium distribution of the spins is given by P{S\J) = 
■exp{-PsH{S\J)), where /3s = Tg \ the Hamil- 

^ J2i<i JijSiSj, and 



Zs{J 

Ionian is given by H{S\J ; ~ ^i^j 
Zs{J) is a partition function under a given J. Within 
the long evolutionary time scale for J, the spin con- 
figuration driven by a Hamiltonian Hj reaches thermal 
equilibrium. The distribution function of J is given by 
P[J) = ^ exp(-/3jiJ,7(J)), where (3j = TJ^ and Zj is 
the total partition function. The function Hj is generally 
expressed in terms of equilibrium quantities of S and a 
bare distribution Pq{J). Here we set the Hamiltonian of 



2 



J as 



Hj{J)^-^{J)-TjlogPo{J), 



(1) 



where ^'(J) is a fitness function. Tlie bare distribution 
is given by Po(J) = n.<j 7^ ("4/(2^0 )) , witli 
a unit Jo of the interaction. We assume that fitness is 
determined by a specific configuration of given t spins, 
called target spins here. (For example, protein function 
depends on the conformation of a set of residues, and is 
indeed modeled by the configurations of target spins in 
p^). More specificically, we assume that a functional 
phenotype is generated when the configurations of target 
spins satisfy ^* Si ~ tfi with fi being a constant value. 
The remaining N — t spins, called non-target spins, have 
no direct influence on the selection of individuals. The 
fitness function is thus defined by 

vI/( J) = log (m, - J2 ^0 ) ^ log(^'(^))' (2) 
1=1 

where 6 is Kronecker's delta and (• • • ) is the thermal av- 
erage with respect to the spin variables according to the 
equilibrium distribution. The fitness function "^( J) im- 
plies a logarithmic probability for the magnetization of 
t-spins to take the value fj, in equilibrium. Note that it 
does not matter which spins are chosen as targets because 
the model is a fully-connected mean-field model. The 
configuration of t-spins is not important either, because 
of the gauge symmetry, which guarantees that a system 
with any configuration of t-spins can be transformed into 
the system studied here, without altering the thermody- 
namic properties (l7j . The equilibrium distribution of J 
and the total partition function are written as 



p{j)^±Po{jms)f-', 

^,7 



zj = m^-']o, (3) 



where [■■ ■]o means the average over J with respect to 
the bare distribution Pa{J). When Tj = 00 or i = 0, 
the distribution P{J) is identical to Po{J) irrespective 
of their fitness values, and the system corresponds to the 
Sherrington-Kirkpatrick (SK) model. For finite Tj and 
i, the interactions J that frequently lead to the spin con- 
figuration with X)i=i ~ appear with higher proba- 
bility. In this sense, the temperature Tj plays the role of 
the selection pressure in genotypic evolution. 

Assuming that /3j is a positive integer, the quantity 
{'ip{S))^-' can be expressed in terms of f3j real replicas. 
Following the replica method [ll], the total partition 
function Zj can be expressed as 



„ P.! 

Zj^ lim / DJPoi J) Tr TT ^ae"^^ ^"=1 ^■ 



(4) 



where = '0(5'") and = H{S°'\J). The right hand 
side of eq. ^ is originally calculated for a positive integer 
n and /?,/ while keeping /3j smaller than n, and then the 



partition function Zj is analytically continued to non- 
integer Pj and non- integer n with the limit to 0. After 
some calculations, the total free energy can be derived 
as a function of replicated order parameters {qa/3}, their 
conjugate parameters {qq^}, and parameters {jj-a} con- 
jugated with fi, which are determined by self-consistent 
equations. The replicas from the first to /3j-th are sub- 
jected to the external field ipa, and the others are not. 
Taking the difference in the replicas into account, we in- 
troduce a replica symmetric (RS) assumption for as 

qi, if «< A/,/3< A/ 
Qap^ { 92, if a < /?,/, /3 > /?,/ or a > f3j, (3 < f3j (5) 
53, if a> Pj,P> Pj. 

For the conjugate parameters {A"}, it is assumed that 
pLa = pi for any a < /?,/. With these assumptions, the RS 
total free energy density /rs is given by 



/rs(Ts, Tj,^i)= -pn/3jfj. - 
Pi fl3j{l3j~l)qf 



2 V 2 

+ {l-p) logS(0)+plogS(/l), 



(6) 



where p = t/ N . Here S(w) is defined as a normalization 
constant of the distribution 



P{u, V] w) 



' cosh(w -I- y/qiu) \ 



2'kE{w) V coshW{u,v) 



(7) 



where W{u, v) 



91 — V 91, 92, and qs 

arc the conjugate parameters of qi, 92, and q^, respec- 
tively. At fi = 0, the free energy is identical to that of 
the SK model under the RS ansatz. The self-consistent 
equations for the order parameters qi, (72, and 53 are 
given by 

qi = (1 -p)(tanh^(V^u))o 

-hp(tanh^(/i+ Vgiw)>A (8) 
(?2 = (1 — p)(tanh(-\/^M) tanhM^(M, v))o 

+ p(tanh(/i + ■\/^")tanhW^(ti,v))^ (9) 

93 = (1 -p)(tanh^iy(M,u))o +p(tanh^M^(M,u))/i, (10) 

where {■ ■ ■)x denotes the average according to the distri- 
bution ^ at w = X. The conjugate parameters of qiS 
are given by qi ~ P^qi {i ~ 1, 2, 3). The first and second 
terms of the order parameters come from the non-target 
spins and the target spins, respectively. Thus, eqs.(8)- 
(10) can be rewritten as the summation of the non-target 
and the target parts, qi = (1 — p)qf^ + Pql (* = 1,2,3). 
The conjugate parameter p is implicitly determined by 
the equation 



/X = (tanh(/i + \/qiu))p,, 
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where is a given parameter in the fitness function and 
the right hand side depends on jl. The stabihty anal- 
ysis for the RS solutions presented by de Almeida and 
Thouless (AT) [Hi affords three conditions 

ATi = 1-/31(1 -2gi+rn)>0 (12) 
AT2 = {l - /?^(l - (/?j + 4)93 + (/?,/ + 3)r33) } 

X + 1 - Pl({f3j + 1)(1 - 93) + (A/ - l)(gi - ^^22))" 
+ 2/?j(/3j + 2)/3|(g2 " r23)' > (13) 
AT3 = 1 - /3|(1 - 2^3 + r33) > 0, (14) 

where 

rii = (1 -p)(tanh'*(Vli"w))o + p(tanh'*(/i + v/^u))/i 
?'22 = (1 — j5)(tanli^(-\/^u) tanh^ VF(m, u))o 

+ p(tanh^(/i + V^i") tanh^ Vt^(u, v))^ 
'^23 = (1 — _p)(tanh(-\/^M) tanh^ M^(w, u))o 

+ p(tanh(/i + v^"") tanh'^ Vt^(u, w))^ 
^33 = (1 - P)(tanh4 W(u, w))o + p(tanh4 W^(u, w))^. 

We introduce an expectation value for the target mag- 
netization mt — [(X)*=i 'S'i)/^]/3j- When mt — 0, the 
fitness function is also equal to 0. Hence, the adapta- 
tion phase is the region satisfying toj > 0. Following the 
replica method, the target magnetization is given by 



mt = (tanhH/(u,w))^, 



(15) 



which indicates that when 52 = 0, the target magnetiza- 
tion is also 0. Thus, the parameter region with (72 > 
(92 = 0) corresponds to the adaptation (non-adaptation) 
phase, respectively. 

The phase diagram on the T5 — Tj plane at p = 0.2 is 
shown in Fig. [TJ Here we focus our attention on the case 
with /i = 1, and we set /t to be sufficiently large to satisfy 
the self-consistent equation eq. ([TT|) with /i = 1. We de- 
fine the transition temperatures Tj^ and such that 52 
and gs are positive or zero, respectively, while qi takes a 
non-zero value at any finite T5. At Tj > 1, the transition 
temperature Tj'* is equal to 1 and the temperature Tj^ 
is smaller than Tj^ = 1. Adaptation occurs at T5 < Tj^, 
but the AT stability conditions AT2 and AT3 are already 
violated at T5 = 1. A preliminary Monte Carlo simula- 
tion indicates that the transition for (72 > and replica 
symmetry breaking (RSB) occurs at T < Tj'' = 1 [lO]. 
At Tj < 1, coincides with r|\ while RSB occurs at 
a lower temperature at which AT^ — 0. Thus the adap- 
tation phase Ts < consists of RS and RSB phases, 
separated by the line AT3 = 0. The RS adaptation phase 
is thermodynamically stable at Tf^^ < T5 < T^^, where 
T^^^, given by AT3 = 0, is the boundary between the 
RSB and RS phases, and T^^ is the transition tempera- 
ture for 52 and 93, Tg ^ = Tj^ = Tj^. As p decreases, the 
region of the RS adaptation phase becomes narrower and 
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FIG. 1. (color online) Phase diagram on the Ts — Tj plane 
at p = 0.2. □ and O indicate the transition temperatures 
and T|^, respectively. The adaptation phase appears in 
the temperature region lower than □. A and v nia,rk the 
boundary of the RSB phase and RS phase. 



eventually the lines of AT2 0, AT3 = 0, Tj" , and 
merge to T5 = 1 for any Tj aX p = 0. In this limit, the 
present model is identical to the SK model whose spin- 
glass transition with RSB occurs at T5 = 1 independent 
of Tj. 

To distinguish the interactions evolved in the RS phase 
from those in the RSB phase, we calculate the equi- 
librium frustration parameters. Indeed, the frustra- 
tion characterizes the interactions in spin glasses. It is 
defined as the product of JijS along a minimal loop. 
When the interactions among the three spins satisfy 
JijJjkJki < 0, the energy per spin cannot reach the 
minimum valuejand such interactions are said to be 
frustrated [13, [HI- In the present model, the tar- 
get spins play a distinct role because their configura- 
tion determines the fitness function. Hence, by distin- 
guishing the target spins from others, we introduce the 
frustration parameters as $1 = X]i<j<t I'^ul/S/ 

*2 = diN-t) T,i<j<t T,k=t+l[J^kJjk]fiJ, where C| is the 
number of interactions between the target spins. When 
$1 = 0. the interactions between the target spins are ran- 
domly distributed; however, when $1 > 0, ferromagnetic 
interactions are dominant. The ferromagnetic interac- 
tions between the target spins energetically favor the spin 
configuration with TOt = 1. The frustration parameter $2 
is the average correlation of the interactions between the 
target and non-target spins. When the interactions that 
couple a non-target spin Sk to the target spins Si and 
Sj satisfy the condition JikJjk > 0, the target configura- 
tion Si = Sj ~ 1 is stable irrespective of Sk- Therefore, 
the finite frustration parameter $2 > implies that the 
configuration with mt > is energetically supported by 
the interactions between target and non-target interac- 
tions. Under the RS ansatz, the frustration parameters 
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FIG. 2. (color online) Tg dependence of the frustration pa- 
rameters (a) "l"! and (b) $2 at Tj — 0.2. The vertical axis in 
(a) and (b) are rescaled with y/N and A*", respectively. The 
RSB transition temperature, which weakly depends on the 
value of p shown here, is indicated by the dashed lines. 



(16) 



(17) 



Here the coefficients N~-^/^ and N~-^ reflect the change 
in the order of the interactions into 0{N~^^-^) through 
the evolution [22| • As seen in Fig. [21 for any p, the frus- 
tration parameters $1 and $2 increase with a decrease 
in Ts down to T5 = T^^, but with further decrease in 
Ts < Tg ^, $1 and $2 decrease. Thus, the frustration 
is minimal at around the transition temperature T^^. 
This result is consistent with the behavior of the energy 
[20j . The configurations of the interactions that evolved 
in the intermediate temperature range Tg^^ < Ts < Tj^^ 
have smaller frustration in the interactions between tar- 
get spins and those between target and non-target spins. 

In summary, we employed a spin-glass model of adap- 
tive evolution to discuss evolutionary robustness in terms 
of statistical physics. Our analysis showed the existance 



of two kinds of adaptation phases, an RS adaptation 
phase at T^^^ < Ts < ^ and an RSB adaptation 
phase at Ts < Tg-^. The equilibrium properties of the 
interactions were characterized by the frustration param- 
eters, which showed that the RS adaptation phase ener- 
getically supports the target configurations by suppress- 
ing the frustration in the evolved interactions. 

Now we discuss the biological relevance of our results. 
An evolved system in the RS phase is robust to noise 
in the dynamic processes and to genetic change. The 
relaxation dynamics of spins progresses smoothly with- 
out becoming stuck at any metastable states. In the RS 
phase, the adapted phenotype, that is, the target spin 
configuration, is a unique stable state that is reachable 
from any initial conditions after a short time of relax- 
ation. This dynamical process agrees well with that of 
the funnel landscape in protein folding [lTI. [lp, a s is also 
observed in evolution dynamics in biology [l4[l5|. Next, 
the self-averaging property in the RS phase guarantees an 
identical equilibrium distribution of the phenotype even 
if the genotype J is distributed around the evolved point. 
An identical phenotype is generated irrespective of geno- 
typic variance, which is known as genetic canalization 
[23|. However, phenotypic robustness is lost at lower tem- 
peratures by RSB, as represented by the appearance of a 
continuous overlap function. Thus, our findings provide 
an evolutionary interpretation for RSB and also confirm 
a positive role for thermal noise in shaping the funnel-like 
dynamics and robustness to mutation. 

Finally, despite the use of a simple statistical-physics 
model of interacting spins, we expect our findings to hold 
true in other problems involving evolutionary and devel- 
opmental dynamics. In addition, the proposed replica 
formalism could function as a theoretical basis to under- 
stand the evolution of robustness in general. 
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